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Particles dispersed in a liquid crystal above the nematic-isotropic phase transition are wetted 
by a surface-induced corona of paranematic order. Such coronas give rise to pronounced two- 
particle interactions. In this article, we report details on the analytical and numerical study of 
these interactions published recently [Phys. Rev. Lett. 86, 3915 (2001)]. We especially demonstrate 
how for large particle separations the asymptotic form of a Yukawa potential arises. We show that 
the Yukawa potential is a surprisingly good description for the two-particle interactions down to 
distances of the order of the nematic coherence length. Based on this fact, we extend earlier studies 
on a temperature induced flocculation transition in electrostatically stabilized colloidal dispersions 
[Phys. Rev. E 61, 2831 (2000)]. We employ the Yukawa potential to establish a flocculation diagram 
for a much larger range of the electrostatic parameters, namely the surface charge density and the 
Debye screening length. As a new feature, a kinetically stabilized dispersion close to the nematic- 
isotropic phase transition is found. 

PACS numbers: 82.70.Dd, 61.30.Cz, 61.30.Hn 



I. INTRODUCTION 

Colloidal dispersions are suspensions in a host fluid of 
solid or liquid particles, of radius ranging from 10 nm 
to 10 /im [l|. They are long-lived metastable states of 
matter that present great interest both from a fundamen- 
tal point of view and in applications, as, e.g., in paints, 
coatings, foods and drugs. Their stability against floccu- 
lation is a key issue in colloidal physics, since the prop- 
erties of a colloidal dispersion drastically change when a 
transition from a dispersed to an aggregate state occurs: 
to prevent coagulation due to attractive van der Waals 
forces, colloidal particles are usually treated in order to 
produce electrostatic or steric repulsive interactions. 

Recently, great attention was paid to liquid crystal 
colloidal suspensions and emulsions, i.e., dispersions of 
solid particles or liquid droplets — respectively — in a liq- 
uid crystal 12, 121 Id, Ifil • In the case of emulsions in 
a nematic phase, the radial anchoring of the nematic 
molecules at the surface of the droplets yields topolog- 
ical defects of the nematic texture in the vicinity of the 
droplets. These defects produce strong repulsions that 
stabilize the liquid droplets against coalescence 4] . Simi- 
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lar effects are present also in other anisotropic fluid hosts, 
such as lyotropic solutions of anisotropic micelles 0, 
and in different liquid crystal phases, as in cholester- 
ics 0. Such systems form composite materials with un- 
usual properties [^ Il0lll]| . 

Even more recently, there has been a growing interest 
in the pretransitional surface-induced interactions me- 
diated by the paranematic order arising in the vicinity 
of surfaces at temperatures above the nematic-isotropic 
phase transition. In |l2j . the force between two par- 
allel plates immersed in a liquid crystal slightly above 
the bulk isotropic-nematic transition was theoretically 
investigated. In j^, the force between a flat surface 
and a glass microsphere, mediated by the surface-induced 
paranematic order, was experimentally measured using 
an atomic force microscope. The measurements were in- 
terpreted in terms of the force between two flat surfaces 
using the Derjaguin approximation jl4| . 

Interesting effects regarding the stabilization of colloids 
are thus to be expected in the isotropic phase of a ne- 
matogenic material when the colloids are wetted by a 
corona of paranematic order |0, 0, 0|. Two effects 
compete: an attraction due to the favorable overlap- 
ping of the paranematic halos (which reduces the vol- 
ume of the thcrmodynamically unvaforable paranematic 
phase) and a repulsion due to the distortion of the di- 
rector field. The vicinity of a phase transition may give 
a critical character to the stabilization mechanism and 
yield rich reversible phase-separation behaviors, as pre- 
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dieted in 0,^3 ^ simpler system with a scalar order- 
parameter. In Ref. p^. the interaction of two spherical 
particles immersed in a liquid crystal above the nematic- 
isotropic phase transition was investigated using a uni- 
axial ansatz function for the paranematic order, in the 
limit of sphere radii large with respect to the nematic 
coherence length. This approximate solution was used in 
Ref. pO] to analyze the stability of a suspension of such 
colloidal particles, in the presence of destabilizing van der 
Waals attractive interactions and stabilizing electrostatic 
repulsions. In Ref. |l6l |. an exact numerical solution for 
the above paranematic interaction was obtained, using a 
multipolar expansion of the tensor order-parameter and 
taking into account biaxiality. The possibility of stabi- 
lizing colloidal particles was discussed, emphasizing the 
case of rather small particles, i.e., of size comparable to 
several nematic coherence lengths. A large-separation 
asymptotic analytic form of Yukawa type was also ob- 
tained, that was shown to very closely describe the ex- 
act interaction between large particles up to separations 
comparable to the nematic coherence length Alter- 
native derivations of the Yukawa potential were given in 
Ref. ^3 based on a Debye-Hiickel type approximation 
and a geometric view. 

In this work, we revise in detail the solution found in 
Ref. 16], which is exact in the limit of weak surface- 
induced paranematic order. We give a particular empha- 
sis on the meaning and limits of validity of the approxi- 
mate asymptotic solution. The latter is used to analyze 
in detail the stability of a colloidal suspension in the pres- 
ence of van der Waals and electrostatic interactions. We 
obtain and discuss general diagrams for the stability of 
the suspension as a function of the various relevant pa- 
rameters. 

The detailed plan of our paper is the following. In 
Section we describe our Landau-de Gennes model 
for the bulk and the surface free-energy. In Section UTTI 
we present the corresponding equilibrium equations and 
formulate the general multipolar expansion for the bulk 
order-parameter in spherical coordinates. This expansion 
is used in Section Hvl to obtain the exact solution for an 
isolated spherical particle. The interaction of two such 
particles is discussed in Section [V] In particular, in Sec- 
tion IV B II we derive the asymptotic interaction energy. 
The asymptotic solution is analyzed in Section FV B 21 in 
terms of the superposition of single-particle solutions. In 
Section rV B 31 we discuss the limits of validity of the Der- 
jaguin approximation. In Section IV CI we present an ex- 
act numerical solution for the interaction of two particles. 
The comparison between the asymptotic interaction en- 
ergy and the exact numerical one is performed in Sec- 
tion ^^^0 In Section FV C 21 we discuss the texture be- 
tween two interacting particles, and in Section FV C 31 we 
analyze the defect ring that appears in the paranematic 
texture between the two particles. Finally, in Section FVll 
using our asymptotic solution, we discuss the stability of 
such colloidal dispersions. 



II. MODEL 

Let us first describe the nematogenic phase in which 
the colloids will be placed. Nematic liquid crystals are 
anisotropic liquid phases, in which elongated molecules 
display a long-range orientational order. This order 
is described by a symmetric traceless tensorial order- 
parameter Qij = 1,2,3), since nematics are non- 
polar • The eigenvectors of Qij represent the axes of 
main molecular orientation and its eigenvalues describe 
the amount of orientational ordering in each direction. 
Usually, nematics are uniaxial phases, however biaxiality 
naturally arises in inhomogcncous situations, e.g., in the 
vicinity of defects 23] . The Landau-de Gennes expansion 
of the bulk free-energy density has the form |3 125| 

fb — ^ QijQij 2 QijQjkQki ^" ^ i^QijQij) 

+^Q ij,kQij,k ■ (1) 

Summation over repeated indices is implied and comma 
indicates derivation. The coefficients are such that a = 
a{T—T*), with a, c, and L positive. The temperature T* 
is the lowest temperature T at which the isotropic phase 
can exist. The presence of the cubic term QijQjkQki in 
the homogeneous part implies that the nematic-isotropic 
transition is first-order. To simplify the description, 
we have introduced only one gradient term in Eq. 
which corresponds to the usual one-constant approxima- 
tion [2^ . The most general expansion of the surface free- 
energy density is |23 

fs — viQijViUj + QijQij + QijQjk'-'iVk 

+ ^Q^jQkll^r'^jl^k^l +0{Q^), (2) 

where the z^^'s (i = 1,2,3) are the components of the 
outward normal u of the surface. Owing to the linear 
term, the surface always locally favor a nonzero order. 

Since we are dealing with the isotropic phase, and, in 
addition, we assume a weak surface-induced order and 
temperatures not too close to the phase transition, we 
can neglect the third and higher-order terms in ft and 
fs ^,1^. Then, the free-energy being quadratic, exact 
calculations arc feasible. To simplify, we shall also retain 
only the simplest quadratic surface term, by setting V22 — 
and V23 — 0. Then the free-energy densities can be put 
in the form 

fb — "^QijQij ~t" ~^Qij,kQij,ki ('^^) 

/. = \w (q., ~ Qf) {Q^, ~ Qf^) , (3b) 

where 

gg) = S, [v.v, - ^-5,^ (4) 
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is the preferred order-parameter at the surface. The 
surface free-energy density (j3b|l is compatible with 
the experimentally measured anchoring in the nematic 
phase |23 |. 

In the following, unless otherwise specified, we shall 
normalize lengths with respect to the nematic correla- 
tion length ^ = {L/ay/^ and energies with respect to 
Fq = aS^^Sg. Because of the conditions that Qij must be 
symmetric and traceless, the free-energy densities Q are 
not diagonal in the five independent components of Qij . 
Defining the following set of components: 



(5a) 
(5b) 
(5c) 
(5d) 
(5e) 
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the normalized free-energy takes the following diagonal 
form: 



F 



4 
■i=0 



9- + {^q^r 



- w 



(9.-9^)^ 



(6) 

where cq = 1/4, ci = C2 = C3 — 1, C4 = 3/4, and 
w = W/a^ is the normalized anchoring strength. The 
first integral is over the bulk and the second one over 
the surface of each colloidal particle. Thanks to this 
transformation, the tensorial problem associated with the 
paranematic order between the colloids is reduced to the 
superposition of five independent scalar problems. 



III. EQUILIBRIUM EQUATIONS AND 
SOLUTIONS 



operator by — 1. We shall therefore use multipolar 
expansions. Actually, the operator — 1 is the opera- 
tor associated with a massive boson field, which implies 
a short-range interaction [2^ . 



A. Multipolar expansion 

In spherical coordinates (r, 9, 0), the spherical harmon- 
ics Yejn{0,(f>) (for their definition see Appendix are 
eigenfunctions of the angular part of the Laplacian op- 
erator and form a complete basis. Then, each of the 
Qi^s can be expanded as 

00 e 

q,{r, 9,ct>)^J2 E Mr)Yi„,iO, 0), (9) 

e=o m=-l 

where the functions ui{r) obey the equations 

^(.^^) = p + ,(,+ l)]„„ ,10) 

which follows by inserting expansion (jSJ into the bulk 
equation (0 . The solution of Eq. (|10|l that is regular at 
infinity can be expressed as 



(11) 



where the Kgj^xj^k''^ are modified Bessel functions of half- 
integer order. Explicitly, 



uo{r) = 



r 



Mr) - — 1 + - , 
r \ r J 

e-"^ / 3 3 
"2W = — ( l + - + ;3 



(12a) 
(12b) 
(12c) 



21+1 

ue+i{r) = ue(r) + ue-i{r) for £>!. (12d) 



At equilibrium, the order-parameter texture minimizes 
the free-energy By setting to zero the variation SF 
of the free-energy, associated with arbitrary infinitesimal 
variations of the qi components, we obtain the equilib- 
rium equations 



{V - l)q, = 0, 



in the bulk, and 



(9.-9^) 



(7) 



(8) 



on the surface of each colloidal particle. Determining the 
paranematic order outside the colloidal particles is there- 
fore somewhat similar to solving an electrostatic problem 
with mixed boundary conditions in which the potential 
is replaced by the qi fields and the standard Laplacian 



With such an expansion, the paranematic order is thus 
completely described by the set of coefhcients gf'". 



B. Equilibrium free-energy 

By integrating by parts Eq. ^ and using the equi- 
librium equations Q) and ©, the normalized free- 
energy (jSJ associated with an equilibrium solution can 
be written as 

f = «;^cJ9f (9f -9.)dV, (13) 

where the integrations run over the surfaces of all the 
particles. 
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IV. PARANEMATIC ORDER AROUND ONE 
PARTICLE 

Let us consider a spherical particle of reduced ra- 
dius R. In spherical coordinates {r^9,(j)), the normal to 
the surface of the particle is given by u = sin 9 cos (/> x + 
sin0sin(/)y -f COS0Z. Therefore, the preferred surface 
order-parameter Q has the following g-components [see 
Eqs. ©] 



1.0 



sin 9 cos 20 



(0) 



(0) 
12 



sin 26* sin (/) = 

sin 26* cos <f> ~ 

■ sin^ 9 sin 20 
1 



cos^e-- = -Y2oi9,(b), 




where 



yR 



— 



1 

2 
1 



(15a) 
(15b) 



are the real and imaginary part, respectively, of the 
spherical harmonics. The boundary conditions ^ be- 
come simply 



dr 



w 



r=R 



q,{R,9,q^)-ql°\9,q^) 



(16) 



Replacing in Eq. qi by its multipolar expansion 
and gl"'' by its expression given in Eqs. (|14() . and iden- 
tifying the coefficients of the spherical harmonics, yields 
the solution 
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\ 



2A{R, w)e^U2ir), 



Y2o{eA)J 



where the amplitude A{R, w) is given by 

R'^w 

~ R^{w + 1) + R^{3w + 4) + 3R{w + 3) + 9 



(17) 
(18) 



and U2{r) is given by Eq. H12c|l . 

Since Eq. (|17|l has the same angular structure as 
Eqs. (|14() . the order-parameter is uniaxial with a radial 
director n — i.e., Qij — S{r){viVj — -j^y), where 



S{r) ^ SoA{R,w) 1 + - + — 



3 3 \ e^-" 



(19) 
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FIG. 1: Normalized scalar order-parameter S{R)/So at the 
surface of one isolated spherical colloid as a function of its 
normalized radius R for various reduced anchoring strength w. 



is the uniaxial scalar order-parameter. Note that what- 
ever the size of the particle, the paranematic order re- 
laxes on the coherence length ^ (equal to unity in reduced 
units). The behavior of the scalar order-parameter at the 
surface of the particle {r = R), as a function of the ra- 
dius R, is shown in Fig. ^ For large values of R, the 
surface order tends to S{R) ~ 50^/(1 -I- w) which cor- 
responds to the limit of a planar interface. For i? < 10 
(in units of ^), the surface order-parameter is reduced 
because of the energy cost associated with the splay of 
the director. 

The free-energy of the particle is obtained from 
Eq. lO. With the help of Eqs. lO and (HTjl and the 
orthonormality of the spherical harmonics, we find 



Fi = —A{R,w) 



R 



9 9 



(20) 



V. INTERACTION OF TWO PARTICLES 

Let us now consider two identical spherical particles 
of reduced radius R separated by a center-to-center dis- 
tance d. We introduce three spherical coordinate sys- 
tems: a global one, (r,9,(l)), symmetrically placed with 
respect to the particles, and two local ones, {ri,9i,(j)) and 
(r2,02,0), centered on the two particles, as indicated in 
Fig. El 



A. General formalism 

Using the symmetry of the system allows to reduce 
the number of free parameters in the tensorial order- 
parameter field. For a given point P, the plane Hp pass- 
ing through P and the z-axis is a symmetry plane. Thus, 
the frame (e^'^^ e(^\ e^^^) in which Qij is diagonal has two 
directions e'^' and e^^^ in the plane Hp and the third, 
e'-'^^ normal to it. This frame can be parametrized by 
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t Z2 



FIG. 2: Geometry for the calculation of the paranematic in- 
teraction between two identical spherical colloidal particles 
separated by a center-to-center distance d. The origin of the 
z-axis is in the midpoint O between the centers of the two 
particles, that of z\ in the center Oi of the upper particle, 
and that of 22 in the center O2 of the lower one. 



the angle Q{r,&) that e*^^) makes with the z-axis: 

e*-^-* = sin8cos0x + sin8sin0y + cos9z, (21a) 



3(2) ^ 

.(3) _ 



COS COS 0x + COS sin^y — sin0 z, (21b) 

(21c) 



sm (px -f cos (py. 



Therefore, the order-parameter tensor Qij can be written 
as 

Q,^=Aie«e«+A.efef +A3ef)ef, (22) 
where its eigenvalues satisfy 

A3 = — (Ai + A2), 



(23) 



since Qij is traceless. The corresponding q- 
components (O are therefore 



go = a(r, 6*) cos2(^, 
qx = /3(r, 6*) sin 0, 
qi = I3{r,0) cos 4>, 
1 
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-a{r, 9) sin 20, 



(24a) 
(24b) 
(24c) 

(24d) 

(24e) 



with 



7(r,0) 



^ (1 + sin^ e) + ^ (1 + cos^ e) , (25a) 

Do Do 



sin 26, 



Ai — A2 
2^0 

Aicos^e + ^sin^e. 

Do DO 



(25b) 
(25c) 



Hence the order-parameter field around two particles is 
fully described by the three fields a{r,6), P{r,6), and 
lir,0). 



1. Multipolar development 

For two particles, the most general expression of Qij 
in terms of multipoles must be a sum of two multipolar 
expansions of the kind each one centered on each 
one of the particles. The reason is the following. Since a 
multipolar expansion such as (|5J) is singular at the origin 
(r = 0), the latter must be put inside one of the par- 
ticles. Moreover, using only one expansion, one would 
impose that the analytical continuation of the texture is 
regular inside the other particle (because the multipolar 
expansion is regular everywhere but at r = 0), which is 
not required. Therefore, in the most general case, there 
must be a second multipolar source inside the other par- 
ticle. For two identical particles, this is also obvious from 
symmetry considerations. 

Taking into account the general form (|24|l and the fact 
that (x, y) is a symmetry plane, the double multipolar 
expansion takes the form (see also Appendix^: 



a{r,9) = ^^^^afUfXrp)Pi>^{cos9p), (26a) 

1=2 p=l 
00 2 

/3(r,0) = ^^(-l)P/3,?/,(rp)P/(cos0p), (26b) 
1=1 p=i 
00 2 

7(r,0) = ^^7,u,(rp)P/(cos0p). (26c) 

1=0 p=i 

Thus, at this point, the whole order-parameter texture 
is uniquely determined by the set of coefficients ae (with 
i = 2,3,...), Pi (with e = 1,2,...), and (with £ = 
0,1,...). 



2. Boundary conditions 

Because of the presence of the (x, y) symmetry plane, 
we need to impose the boundary conditions ^ only on 
one of the particles, let say particle 1. These boundary 
conditions on the g-components transform to boundary 
conditions on the fields a, (3 and 7 [see Eqs. (I24|) ]. which 
in the coordinates (ri, Oi) relative to particle 1, take the 
form 



da 
dri 
d(3 
dri 

dri 



ri=R 



ri=R 



ri=R 



«(ri=i?,0i)-ip2'(cos0i) 
/3(ri=i?,0i) + ip2'(cos0i) 
7(ri=P,^i)-^P2"(cos0i) 



, (27a) 
,(27b) 
(27c) 



Note that these equations are decoupled; therefore the 
three fields a, (3 and 7 can be treated independently. 
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3. Coordinate transformation 

To impose the boundary conditions (|27|l to the muhi- 
polar expansions H2t)|l . we need to express the spherical 
coordinates r2 and 62 in terms of ri and 61. From Fig. [21 
we easily obtain 



r2 = V n 
cos 62 = 



^(P + 2dri cos6'i, 
ri cos di + d 



+ 2dri cos 6*1 



(28a) 
(28b) 



mined analytically, owing to the intricate relations H28|l 
that link the local coordinate systems centered on the 
two particles. However, in the limit 0? 3> 1, i.e., for 
distances large with respect to the nematic coherence 
length ^ (even if c? — 2R <C i?) , we can obtain an asymp- 
totic solution by expanding the unknown multipolar co- 
efficients a£, /S^, and 7^ in series of e~'^/c?" up to a given 
order n. In the following, we shall obtain the lowest or- 
der expansion, n — \. It turns out that to this order the 
only non-zero terms in the multipolar expansions l|26(l are 
those up to ^ = 2. 



B. Asymptotic solution 

The exact two-particle order-parameter texture H26|) 
satisfying the boundary conditions (|27|l cannot be deter- 



To begin with, we illustrate the strategy for comput- 
ing the 7£ coefficients. We start by writing the expan- 
sion H2t)cll on particle 1 (ri = i?) up to ^ = 2: 



7(ri =i?,0i) 



70 [uo{R)Po\cose^) + u^{R2)Pa\cosQ2)\ + 7i [iii(i?)Pi°(cos0i) ui(i?2)Pi"(cose2)] 

72 [u2{R)P2\cOSe^)+U2{R2)P2°{cOsQ2)\ , 



(29) 



where R2 = R2{0i) is the r2 coordinate H28a|l evalu- 
ated on particle 1 (ri — R), and, similarly, cos 82 = 
cos 82(^1) is the cos 6*2 coordinate Ij28b|) evaluated on 
particle 1. Next, we asymptotically develop the functions 
Mn(^2)^n°(cos 82) up to tcrms of order e^'^/d: 



and we redevelop this expansion on the complete basis 
of Legendre functions of first kind P„°(cos0i), using the 
orthogonality relations (|A3I) . Truncating the expansion 
at ^ = 2, for the sake of consistency, we arrive at 



u„(i?2)P„°(cose2) ~ (_i)»£^e-^™^^\ 

d 



(30) 



e ^ f sinhi? 
sinh R 3 cosh R 3 sinh R 



,(i?2)P„°(cOSe2) ~ (-1)"— |^^^^Po"(cOS0i)-3 



R 



i?2 



i?3 



cosh R sinh R 



R i?2 

P2°(COS01) 



Pi°(cos6'i) 



(31) 



Replacing Eq. (|31|l into Eq. H29|) yields the required expansion of 7. To match the boundary conditions, let us 
determine the corresponding expansion for the radial derivative of 7. To this aim, it is most efficient to take into 
account that the radial functions l|12() obey the relations 



du„{r) _ nunjr) 
dr r 



(32) 



We then obtain 



97 
dri 



ri=R 



- 70 <j -wi(i?)Fo°(cos^i) + — [Mo(r2)Po°(cos02)],^^fl 



71 

72 



ui{R) 
R 

2u2{R) 
R 



-U2{R) 

-usiR) 



Pi "(cos 01 ) + ^ [«i(r2)Pi°(cOS(?2)],^^fl| 
P2°(COS01) + ^ [u2(r2)P2"(cOS02)], _ 



(33) 
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in which, according to Eq. (|31|l . 

d 



dri 



K(r2)P„°(cos02)]^^^^ ^ (-1)' 



cosh i? sinh R 



R 



i?2 



Po"(cos0i) 



sinh R 2 cosh i? 2 sinh i? 



i?2 



R^ 



Pi°(cos6li) 



cosh i? 4 sinh R 9 cosh i? 9 sinh R 



R 



i?2 



i?3 



i?4 



P2°(COS01; 



(34) 



Finally, we insert the relations (|29|l . (|31|l . (|33|l . and H34|l in the boundary condition H27c|l . Because of the orthogonality 
relations HA3|I . we get three linear coupled equations for the unknown coefhcients 70, 71, and 72. Solving them to 
leading order in e~'^ jd yields, after some lengthy algebra, 



70 



71 



72 



e^R^w {i? (u; + 1) + 1 - e^^ \R (w - 1) + 1]} 



i?4 (w + \f + i?3 (3w2 + 8ui + 5) + i?2 (3iy2 _^ i5y^ _^ _^ gi? (2ui + 3) + 9 

J-'' e-^i?*!;; (u; + 1) + 4. 2) + 2 + e^^ - 1) - i? (w - 2) - 2] } 

R^{w + \f + 2i?4(2w2 + 5w + 3) + i?3(6w2 + 24w+19) + i?2 (3^2 + 30^ + 35) + 3i?(5w + 12) + 18' 



\A(R,w)e^ 



2~d 



1 - 



e2«yt(i?, w) (w - 1) - i?2 (3w - 4) + 3i? (w - 3) + 9] 



RH 



(35a) 



(35b) 



(35c) 



The calculation of the remaining coefficients ag and /?£ 
is much simpler since, to leading order in e~'^ /d^ they do 
not depend on d, because of the asymptotic behavior 



Wn(-R2)^'n^(c0Se2) CX 



d2 ' 



Un{R2)Pn (cos 82) OC 



(36a) 
(36b) 



These coefficients are therefore given by their expressions 
for an isolated particle 



a2 ~ ^A{R,w)e^, 
Pi ^ 0, 

p2 ^ -i^(i?,«;)e^. 



1. Asymptotic interaction energy 



(37a) 
(37b) 
(37c) 



To calculate the free-energy we use Eq. ((T^ . By sym- 
metry, it is equal to twice the contribution corresponding 
to the surface of particle 1. In order to calculate this con- 
tribution, we express a{ri = R,Oi), /?(ri = R,Oi) and 
7(ri — R,9i) as an expansion in Legendre functions: 

00 

a(ri =i?,6'i) = ^A„P„2(cos6»i), (38a) 

n=2 

00 

/3(ri =i?,0i) = ^B„P„i(cos0i), (38b) 

n=l 

00 

7(ri =i?,0i) = ^C„P„"(cos0i). (38c) 



Using Eq. p3|) and the orthogonality of the Legendre 
functions, we obtain 



F = —ttR^w (10 - 1274, + I2B2 - 3C2) . 
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(39) 



The coefficients A2 , B2 and C2 are easily calculated from 
the asymptotic expansions of a, (3, and 7 given in the 
preceding section. We arrive at P = 2 i^i 
Fi is the free-energy (|20ll of an isolated particle and 



■FiS: where 



poo _ 



-8^i?8w2e-('*-2«) 



M [i?3(w+l) + R^{2>w+A) + 'iR{w+?,)+9\ 



8^ ,2 e"('^-2«) 
' 3 d 



(40) 



n=0 



is the asymptotic interaction free-energy. It is always 
attractive and has the form of a Yukawa potential {29l |. 
In Subsection IV C II using a numerical calculation, we 
shall show that this asymptotic expression is quite good 
even up to separations of the order of ^ or less. 



2. Analysis of the effective Yukawa potential 

There exists an instructive geometric view on the ori- 
gin of the attractive part of the two-particle interaction 
mediated by the surface-induced paranematic order. It 
reproduces the form of the Yukawa potential of Eq. H40|l 
under the assumption of large particles {R^ 1) and was 
presented by one of the authors in Ref. The basic 

idea is that the overlapping of the paranematic coronas 
of the two particles reduces the volume of the energeti- 
cally unfavorable liquid crystal ordering, as illustrated in 
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FIG. 3: The attractive paranematic interaction between the 
particles comes from the overlapping of the two paranematic 
coronas of thickness ^. In the geometric view, the interaction 
energy is equal to the negative of the free-energy of the excess 
orientational order in the dark shaded region [for an exact 
definition see Eq. 14111 1 . The half spaces Vi and V2 are defined 
by the midplane of particle 1 and 2. 



Fig. 13 and therefore the total free-energy. The interac- 
tion energy is the free-energy of the removed orientational 
order. If we denote by the single-particle order- 

parameter field centered on particle n and switch to the 
representation by the components g^-"'' (« = 0, . . . , 4) [see 
Eq. H17|) ]. the interaction energy is calculated as 



int 



V2 



d'7 



(41) 



where V2 denotes the half space of particle 2 (see Fig.|2J). 
This definition leads to the Yukawa potential of Eq. H4U|) 
in the limit R ^ 1, however with half the strength. 
Its advantage lies in its semiquantitative agreement with 
Eq. H4l)|) and its simplicity. Using a type of Voronoi cell 
construction, it is extensible to multibody interactions 
which are important in the study of particle aggregation 
and the formation of ordered crystalline structures. 

It was also demonstrated in Ref. that the Yukawa 
potential of Eq. (jlUl) can exactly be derived by approx- 
imating the tensorial order-parameter through a linear 
superposition of the two single-particle solutions centered 
on particle 1 and 2 



Qvir)^QrAr) + Q^'\r) . 



(42) 



This approach is in full analogy to the treatment of the 
two-particle potential in electrostatically stabilized col- 
loids based on the Debye-Hiickel approximation. 



3. Derjaguin approximation 

The force between two interacting spheres is often cal- 
culated in terms of the interaction energy per unit sur- 
face £{s) between two parallel plates a distance s apart, 
using the so-called Derjaguin approximation jl4j . For 
two equal spheres of radii R a distance s apart, the 
force in the Derjaguin approximation is simply J-d{s) = 



0.6 
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FIG. 4: Relative error D of the force in the Derjaguin approx- 
imation with respect to the exact numerical one, calculated 
according to the procedure outlined in Section lV CI The error 
is plotted as a function of the distance to contact s = d — 2R 
for a particle of reduced radius R = 20 and for various values 
of the reduced anchoring w. Contrary to what expected, the 
approximation does not improve at short distances. 



ttR£{s) [13. Such an approximation is valid if the inter- 
actions are additive and if the radii of the particles are 
large with respect to both the range of the interaction 
and the minimum separation distance s between the par- 
ticles. Although in our case the interaction energy is not 
pairwise-additive, let us examine whether there exists a 
regime in which the Derjaguin approximation holds. 

The interaction energy between two parallel plates, im- 
posing homeotropic boundary conditions, immersed in a 
nematic liquid crystal was calculated in [T^ . For weak 
induced paranematic order in the isotropic phase, we can 
use our quadratic free-energy ©. The resulting order- 
parameter is everywhere uniaxial, with a nematic direc- 
tor everywhere orthogonal to the surfaces. In our reduced 
units, a straightforward calculation gives the exact inter- 
action energy per unit surface 



3(1 -I- u;) [e'' (l + w) + w-iy 



(43) 



that is everywhere attractive. The Derjaguin approxima- 
tion gives the interaction force 



AttRw'^ 



3 (l + w) [e« (1 



w 



■ w - 



(44) 



To check the Derjaguin approximation, we need to 
compare the forces and not the interaction energy. In- 
deed, since the Derjaguin approximation does not hold 
for separations large with respect to the radii of the parti- 
cle, we cannot fix the correct zero level for the interaction 
energy. We calculate the exact force ^ = —dFi^t/dd by 
numerically evaluating the interaction energy F-mt , as will 



9 



10" I 

|A| z 



10 = 



10"^ = 



10" 



= 1 1 1 U 'i ' " 

^^=^=■■==-=--,4^^ 0.1 


1 1 1 Iz 


1 -_7i; = 0.5>^ 


lllllll 


1 1 1 1 


o- — 

\>vvu; = 10 ^ 

\ 

\ \\\ = 
\ \v = 

\ 






w ^1 ^ 


1 


1 \'i \r 



0.1 



1.0 



d-2R 



FIG. 5: Absolute value of the relative error A of the asymp- 
totic interaction energy as a function of the distance to con- 
tact d — 2R for R = 2 and various values of the reduced 
anchoring w. The solid lines indicate that A > 0, the dashed 
lines that A < 0. 



be explained in Sec. IV Cl In Fig. 0] we plot the relative 
error D, defined through — Td (1 + D). As it is seen, 
the Derjaguin approximation is rather crude. Moreover, 
it does not improve as the separation distance s tends 
to zero, as is normally expected. This is most probably 
due to the non pairwise-additive character of the force. 
The error is actually lowest in the intermediate range 
1 <C s ^ i?. Indeed, in the limit s ^ 1 the Derjaguin 
approximation (|44|l becomes 



inRw e ' 
3(H-u;)^ 



(45) 



which coincides, in the limit i? 3> 1, with the exact 
asymptotic force —dF-^^/dd obtained from the asymp- 
totic interaction energy H40() . with s — d — 2R. Care 
should therefore be taken in using the Derjaguin approx- 
imation for distances to contact comparable with the ne- 
matic coherence length ^ [l^. 



C. Numerical results 

When the two colloidal particles are an arbitrary dis- 
tance d apart, we determine numerically the equilibrium 
order-parameter texture. To this aim, we truncate the 
expansions H26|l at some order i = ^max- Then, with 
the aid of the coordinate transformation H28() and of the 
orthogonality relations l|A3(l . we project numerically the 
boundary conditions H27|l onto the Legcndrc functions 
of the first kind P™{9i), with m 2 for the func- 
tion a{r,0), m — 1 for (3{r,9) and m = for 7(r, 6'), 
and £ = m, . . . , ^max- This gives three decoupled sets of 
linear equations that we solve numerically to determine 
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FIG. 6: As Fig.Elbut for spheres of radius R = 20. 



(with £ = l,...,^max), and 7£ (with £ = 0, ...,£max). 
Finally, we vary the order i^max of the expansions to ob- 
tain a given accuracy. Typically, only a few multipoles 
(^niax — 10 20) are needed to reach relative errors of 
the order of 10~^ or less. 



1. Interaction energy 

We determine the exact interaction energy of the col- 
loidal particles by inserting the numerically determined 
expansions 12()l) in the equilibrium free-energy (|13|l . The 
behavior of the exact interaction energy was studied in 
some detail in |1(t| |. Here we concentrate on its compari- 
son with the asymptotic expression (|40|l . to validate the 
latter. We set 



Fint = (1 -f A) , 



(46) 



the unknown coefficients ai (with £ = 2, 



where A measures the relative error from the asymptotic 
solution. In Fig. [S] we show the behavior of A as a func- 
tion of the distance to contact d ~ 2R for small particles 
{R — 2) and various anchoring strength w. As it is seen, 
the asymptotic approximation is quite good up to dis- 
tances to contact of the order of For low anchoring 
strength, w < 1, the actual interaction energy is larger 
than the asymptotic value. For high anchoring strength, 
w > 1, the actual interaction is smaller. The best agree- 
ment is for anchoring strength of the order of w = 1 , for 
which the relative error increases up to 0.6 at particle 
contact. 

For relatively large particles (i? = 20), the relative er- 
ror A becomes even smaller, as shown in Fig. |H| It is 
remarkable that around w = 1 the asymptotic approx- 
imation gives extremely good results up to contact for 
larger spheres. 
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2. Texture around two particles 

According to (|25|1 and H23() , the eigenvalues associated 
with the eigenvectors (|21f) of the order-parameter tensor 

I 



Far from the particles, or for an isolated particle, the 
texture is uniaxial, with A2 — A3. The paranematic 
director n thus coincides with the eigenvector e^-'^^ [see 
Eq. H21a|l ] associated with Ai, which is the eigenvalue 
of largest absolute value. Therefore, by continuity, we 
choose as paranematic director n the direction e'-'^' . In 
general, however, the tensor order-parameter Q is bi- 
axial and is characterized by the uniaxial scalar order- 
parameter S in the direction of the paranematic di- 
rector and by the biaxiality B in the plane orthogo- 
nal to it. They are defined by setting in the diago- 
nal frame (e^^^ e^^^ e^'^^), in which Q = diag(Ai, A2, A3), 
Q = diag{2/3 S, -1/3 S + B, -1/3 S-B) Therefore: 

S{r,9) = ^Ai(r,0), (49a) 

B{r,9) = ^[X,{r,9) + 2X2{r,9)]. (49b) 

In Fig. [7| we have plotted the contour lines of the 
scalar order-parameters S (top) and B (middle), and 
the field-lines of the paranematic director n (bottom) 
for small spheres, i.e., R — 2. In Fig. |Slwe have plot- 
ted the same parameters for large spheres, i.e., R = 20. 
When the two spheres are close to each other, the scalar 
order-parameter makes a diffuse halo around the two 
spheres for the small spheres, whereas it is significantly 
nonzero only in the gap between the spheres for the large 
ones. The biaxiality is maximum in the vicinity of a ring 
which defines a defect, of topological charge —1/2, of the 
paranematic director n. The defect location corresponds 
to the condition S = B, i.e., Ai = A2 [see Eqs. (|^ ]: the 
two orthogonal direction e*^^) and e^^^ in the symmetry 



are 



(47a) 



(47b) 
(47c) 



(48) 

I 

plane Hp (the plane of Fig. |21 containing the ^;-axis and 
the generic point P) become equivalent. Strictly speak- 
ing, at the defect location the texture becomes again 
uniaxial, with a paranematic director orthogonal to the 
Hp-plane and a discotic-like order. The main difference 
between the small and the large spheres is that the de- 
fect is in the region where S is large in the former case, 
whereas it is located in a region where S is almost zero 
in the latter. Regarding the texture, one can see that the 
field lines are rather smoothly bent in the case of small 
spheres, while for large ones they start almost radially 
from the particles and present a kink in the mid-plane, 
in the region where S is small. 



3. Defect position 

The defect ring is located in the plane z = of Fig. |21 
On this plane, /3 = for symmetry reasons. Then, the 
location of the defect, which corresponds to the condition 
Ai = A2, according to Eqs. (|47a|l and H47b|l is such that 
a{r, 9 — tt/2) = 3'^{r, 9 — tt/2). Let us calculate asymp- 
totically, for d ^ 1, the defect position. At zeroth order, 
the total Q tensor is the superposition of the Q tensors 
of the two isolated particles. Then, one easily finds that 
the defect is located at a radial distance h — d/2 from 
the axis joining the centers of the particles. This corre- 
sponds to a defect ring having a diameter equal to the 
distance between the centers of the two spheres. 

To verify this lowest-order prediction, in Fig.|51we plot 
the ratio between the ring diameter 2h and the distance 
between the two spheres as a function of the distance to 
contact d— 2R. As one sees, the approximate estimation 



Ai(r,0) ^ oMHtM) 
^0 4 

A2(r, 9) ^ ajr, 9)+^{r, 9) 

So 4 
X3{r,9) a{r,9)+j{r,9) 



air, 9)-3"fir,9) 



PHr,9), 



a{r, 9) — 37(r, ( 



/52(r,0), 



5*0 2 

while the angle 8(r, 9) that the eigenvector e*^^^ makes with respect to the z-axis is 



e(r,6l) = tan"^ 



a(r, 9) - 37(r, 9) + \l [a(r, 9) - 37(r, 9)Y + 16/32(r, 9) 
4/3(r,0) 
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FIG. 7: Paranematic order-parameter between two spheres 
with reduced radius R = 2 and reduced anchoring w — 4. 
Top: contour hues of the scalar order-parameter S; middle: 
contour lines of the biaxiality parameter B; bottom: field 
lines of the paranematic director n. 




FIG. 8: Same as Fig. □ but for R = 20. 




2 d-2R ^ 

FIG. 9: Ratio between the defect ring diameter 2h and the 
distance between the spheres d as a function of the distance to 
contact d — 2R of the spheres. Full lines: R = 2; dashed lines: 
R = 5. The two upper curves are for a normalized anchoring 
strength w = 0.1, the two lower ones for w — 4. 



VI. FLOCCULATION TRANSITION 



is rather good: for lower anchoring strength, the defect 
diameter is actually slightly larger, while for higher an- 
choring it is smaller. The deviation from the lowest-order 
solution h = d/2 decreases as the spheres move further 
apart or become bigger. In this case, the defect ring lives 
in a region where the order-parameter is small. 



Section FV C II fsee also Ref. and the Yukawa po- 
tential of Eq. 1401) show that for large particles (i? ^ 1), 
the two-particle interaction mediated by the nematic wet- 
ting layer (in the following abbreviated by Ulc) is dom- 
inated by its strong attractive part. The consequences 
of this feature on the stability of a dispersion of colloidal 
particles were investigated in detail by Borstnik, Stark 
and Zumer in Ref. p(J|. The main idea was to use, e.g., 
the screened electrostatic interaction of charged parti- 
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cles to stabilize a colloidal suspension against the attrac- 
tive van der Waals interaction well above the nematic- 
isotropic phase transition. Then by lowering the temper- 
ature towards the transition temperature Tc, which in- 
creases both the strength and range of the liquid-crystal 
mediated interaction, a flocculation transition can be in- 
duced, i.e., the particles start to aggregate. Whereas 
the transition itself has not been observed so far, the at- 
tractive two-particle potential has been demonstrated in 
excellent experiments by Kocevar and Musevic using an 
atomic force microscope [s^, |^ ■ 

In this article, we study again the flocculation transi- 
tion for two reasons. In Ref. IJoj the interaction poten- 
tial Ulc was calculated by modeling the liquid-crystalline 
order with the help of an ansatz function. A weak repul- 
sive barrier in U lc occurred which seems to be an artifact 
in the construction of the ansatz function. Secondly, by 
using the analytic form of the Yukawa potential for C/lCj 
the calculation of the total interaction energy becomes 
much easier. Therefore, in Section FVIBI we shall discuss 
a flocculation diagram in terms of the relevant parame- 
ters of the electrostatic interaction which covers a much 
larger range than in Ref. (23 ■ We will, however, confirm 
the flocculation diagram of Ref. ^3 ■ 



A. Two-particle interactions 

In the following we consider a colloidal dispersion of 
particles subject to van der Waals {Uw), electrostatic 
(Ue), and liquid crystal induced {Ulc) interactions, 
where Uw and Ue are taken according to the DLVO the- 
ory Is^ . In a first approach, we assume that the sur- 
face induced nematic order has no effect on Uw and UE^ 
so that the total two-particle potential can be written as 
Utot = Ulc + Uw + Ue. 

To be concrete, we consider silica particles of radius 
R = 250 nm dispersed in the liquid crystal 4-n-octyl-4'- 
cyanobiphenyl (8CB). The van der Waals interaction be- 
tween identical particles is always attractive. It reads |33 | 



Uw = -^ 



H 
6 



2i?^ 



d2 - 4i?2 



-t- 



2_R2 
d2" 



In 



rf2 



(50) 

where d is the center-of-mass separation of the parti- 
cles and iJ = 1.1 fceT is the Hamaker constant for silica 
particles suspended in the compound 8CB '2^. Strictly 
speaking, via retardation effects, the Hamaker constant 
also depends on the separation d (s^l- However, in our 
problem this dependence is not crucial [23 ■ Note that 
in the following we do not use reduced units: all physical 
quantities keep their dimensions. 

We stabilize the colloidal dispersion by electrostatic re- 
pulsion. The particles possess an electrical surface charge 
density Qs-, and monovalent salt of concentration Up is 
added to the solvent. This results in the two-particle 



potential [34 



Ue 



£2 



(51) 



The range of the interaction is determined by the Debye 
length = [s2kBT/ {87re'^np)Y^^ , where £2 denotes the 
dielectric constant of the solvent and e is the fundamental 
charge. Note that Eq. (|51|l is derived in the Derjaguin 
approximation under the assumption d — 2i?, k"^ <C R. 
For realistic values of and Qs, see Ref. 20J and 31]. 

To include the liquid-crystal mediated interaction, we 
use the Yukawa potential of Eq. H4U|I . The latter, restor- 
ing to all the physical quantities their dimensions, is given 
explicitly by 



ULc^-fueiSoAr- 



(52) 



where the amplitude A must be calculated according 
to Eq. H18|) . with R substituted by i?/^ and with w = 
W/[La(T — T*)]^/2 |;iig reduced surface coupling con- 
stant. We recall also that the nematic coherence length 
is given by ^ = [L/{a{T -T*))]^^^ . In ScctionrvOit was 
shown that for large particle radii the Yukawa potential 
agrees well with the numerical results down to a distance 
to contact of the order of ^. Even for smaller distances 
it gives a good approximation to the interaction energy. 
We use Eq. H52|) to obtain an estimate for the strength 
of the Yukawa potential in the case R^ 



\ULcid^2R)\ 



An 



■ LRSn 



w + 1 



(53) 



The rational function in the reduced surface coupling 
constant w is monotonously growing like w'^ for small w 
and approaching one at large w. Therefore, we find that 
the strength of the Yukawa potential increases, like its 
range f , when the temperature is lowered towards the 
nematic-isotropic phase transition at Tc- For example, 
for the material parameters chosen immediately below, 
we find that the strength increases by a factor of six when 
cooling down from Tc -I- 10 K to Tc- When discussing the 
effect of the interaction Ulc: we shall use the parameters 
of 8CB [a = 0.12 x 10^ erg/^cm^K), L = 1.8 x lO'^ dyn, 
Tc-T* = 1.3 K, and = 314.8 K], which give a co- 
herence length £,{Tc) — 10.74 nm at the phase transition. 
Furthermore, we choose W — 1.25erg/cm2 and 5*0 = 
0.45, in accord with Ref. 20]. Finally, all interactions in 
the next subsection are referred to the thermal energy 
k^T at room temperature, and, as already mentioned, 
the particle radius is 250 nm. 



B. Flocculation diagram 

When discussing the total interaction energy for vary- 
ing parameters, we will find different shapes of the two- 
particle potential which affect the flocculation transition. 
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FIG. 10: Flocculation diagram in terms of the tunable param- 
eters of the electrostatic interaction, i.e., the surface charge 
density qs and the Debye length The four different re- 

gions I-IV characterize how a colloidal dispersion of charged 
particles reacts on Ulc- The points (a)-(e) refer to the dia- 
grams in Figs. Illj and ll2l 



Particles start to aggregate when the interaction poten- 
tial exhibits a minimum Umin < at finite distances. 
If the potential minimum is shallow, i.e., just several 
k^T deep, a particle doublet will break up again, and 
a phase equilibrium of an aggregated and dispersed state 
occurs |34| . The higher potential energy in the dispersed 
state is compensated by the larger entropy. On the other 
hand, when |C/niin| tends to 10 or even higher values, 
strong attraction occurs that leads to a non-equilibrium 
phase with all the particles aggregated. This corresponds 
to the flocculation transition that we aim to investigate. 
As a further feature of our interaction potentials, we will 
also encounter potential barriers which slow down the 
flocculation or even "kinetically stabilize" the dispersion. 

Obviously, the observation of a flocculation transition 
is a matter of time scale on which the experiments take 
place. The theory of aggregation kinetics determines the 
characteristic time r for the formation of a particle dou- 
blet as HIH 



1 i?^ 



I, 



where the factor / is given by 

T = 2B r ^0 exp([/tot/fcBr) 

which contains the ratio 



ds. 



Do 



2s I R for s < i? 
1 for s > i? 



(54) 



(55) 



(56) 



In Eq. (|54|) . is the volume fraction of the particles in 
the solvent, and Do — ksT / {GirriR) is the single particle 
diffusion coefflcient depending on the solvent viscosity rj. 



For particles that stick together when they meet but oth- 
erwise do not interact (/ = 1), t is given by the time an 
independent particle needs to diffuse a distance R times 
the inverse volume fraction 0. Clearly, in a very dilute 
dispersion (0 small), particles hardly come close to each 
other and doublet formation is very rare. The factor /, 
where s = d — 2R denotes the distance to contact, incor- 
porates influences from the two-particle interaction; in 
the presence of a potential barrier, r increases consider- 
ably even for moderate cf>. The ratio Do/D{s) takes into 
account corrections of the diffusion constant due to hy- 
drodynamic interactions. Note that Smin generally means 
s = 0, i.e., the particles stick together at direct contact. 
If, however, for very small particle separations a strong 
repulsion sets in, Smin is taken as the location of the po- 
tential minimum which is responsible for the formation 
of the particle doublet. 

In the following, we assume a volume fraction cj) — 0.1 
so that for a typical shear viscosity of 0.4 P we arrive 
at t/I = 5 s. We say that a dispersion is kinetically 
stabilized by a potential barrier if r > 1 h, which cor- 
responds to / = 720. With the help of the algebraic 
program Maple, we have numerically calculated the fac- 
tor / for our potential J7tot to check when such a kinetic 
stabilization in the flocculation diagram, to be discussed 
below, sets in. In this way, we have determined the tran- 
sition lines between regions I/III and III/IV in Fig. [TUl 
To gain further insight into the factor /, we perform 
a saddle-point approximation. We replace Utot by its 
harmonic approximation around the maximum J/max of 
the potential barrier. Then we can evaluate the integral 
in Eq. (|55|l when we approximate s by Smax otherwise 
and choose Smin — 0. From the magnitude of the curva- 
ture Co = ds'^\s^.^^, we deduce a characteristic 
length ^ = (2fcBr/co)^^^ and make the approximation 
of replacing Smax by ^ which finally gives the estimate 
/ = exp(J7inax/fcBT)/2. From our criterion for kinetic 
stabilization, / = 720, we find ?7max = 7.3 k^T. Surpris- 
ingly, in most cases the exact calculation of the integral 
in Eq. H55(l gave potential barriers with C/max between 7 
and 8 k^T. This demonstrates that the Boltzmann factor 
is the determining quantity in / and therefore in r. 

We will also encounter two-particle interactions where 
the particle doublet settles into a potential minimum 
without traversing a potential barrier. Here the question 
arises how stable the doublet is. Unfortunately, the famil- 
iar Kramers rate "sg] is not applicable, since it involves 
the curvature of the potential maximum over which a 
particle escapes. In our case this potential maximum 
is located at infinity (J7tot = 0)- It has zero curvature, 
which gives an infinite escape time. The theory for our 
problem formulates the escape time as a double integral 
over the particle separation |3^. that cannot be calcu- 
lated straightforwardly. However, we can give an upper 
bound for the escape time [33, defined as the time idiff 
that a particle needs to diffuse a distance R when leaving 
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distance to contact d-2R [nml 



FIG. 11: Total two-particle potential in units of k^T as a function of distance to contact d — 2R. The potentials are shown for 
different parameters of the electrostatic interaction, as indicated in the inset, and at various temperatures relative to Tc- The 
labels (a)-(d) refer to a location in the flocculation diagram of Fig. IIUI 



a potential well of depth Uq 

^""^ = 6i?oexp(-C/o/fcBT)- ^^^^ 

Again, we consider a particle doublet as stable when ^diff 
exceeds 1 h, which for the parameters already introduced 
results in a potential depth of ?7o = 9 k-oT. Compared to 
the free diffusion time R^/{6Do) = 0.5 s, it enhances idiff 
by a factor of 7200. We used this criterion to establish 
the boundary line between region II and III in Fig. ^| 

In the following, we discuss in detail the flocculation 
diagram of Fig. ^1 (calculated for a particle radius R of 
250 nm) . It is plotted in terms of the tunable parameters 
of the electrostatic interaction, i.e., the surface charge 
density Qs and the Debye length k~^. We distinguish 
four regions: in region I the particles are completely ag- 
gregated without Ulc since the electrostatic interaction 
is not sufficient to stabilize the dispersion against the van 
der Waals force. For large surface charge density and 
interaction range, as in region II, the dispersion is com- 
pletely stabilized even if Ulc is turned on. A potential 
minimum at finite distances does not occur, or it is not 
deep enough, according to our discussion above, to trig- 
ger flocculation. Both regions are separated by areas 111 
and IV. In region 111 flocculation is induced by lowering 
the temperature towards Tc- Different types of interac- 
tion potentials occur which are illustrated in Fig. 1111 for 



various locations in the flocculation diagram of Fig.llOlla- 
beled by (a) to (d). In diagram (a), the graph for T ^ Tc 
exhibits the usual "primary minimum" at short distances 
due to the van der Waals interaction which is followed by 
a "primary potential barrier" for increasing separations 
well known in colloids science. The electrostatic repul- 
sion is just strong enough to stabilize the colloidal sus- 
pension kinetically. Lowering temperature reduces the 
barrier and induces flocculation. Moving down the di- 
agonal from left to right in the flocculation diagram, we 
arrive at location (b). In the corresponding Fig. Ilir b). 
the primary minimum at T ^ Tc is no longer visible 
due to the increase in the surface charge density. To in- 
duce flocculation, one has to cool the dispersion closer 
to Tc- Moving further down the diagonal, a shallow po- 
tential minimum at around 26 nm appears at T ^ Tc, 
as illustrated in Fig. Illf c). which leads to a slight phase 
coexistence. However, a clear non-reversible aggregation 
of the particles will be induced when the temperature is 
decreased. A new feature occurs in diagram lllf d). which 
is located close to the transition line to the stabilized dis- 
persion. Whereas in the previous cases the flocculation 
transition sets in gradually with decreasing temperature, 
here a sudden transition is observable within a few tenths 
of Kelvin. Crossing the transition line from location (d) , 
the minimum at T = Tc becomes more shallow and ulti- 
mately vanishes. Region IV in the flocculation diagram 



15 




10 20 30 40 50 60 70 80 
distance to contact d-2R [nml 



FIG. 12: Total two-particle potential in units of k^T as a 
function of distance to contact d — 2R. The potential is shown 
at various temperatures relative to Tc- The parameters of the 
electrostatic interaction is shown in the inset. The label (e) 
refers to a location in the flocculation diagram of Fig. IIUI 

of Fig.^lidentifies a kinetically stabilized dispersion. As 
illustrated by Fig. 1121 at T = Tc a minimum in the to- 
tal interaction energy is separated from large distances 
by a potential barrier which prevents the formation of 
stable particle doublets. The interaction becomes totally 
repulsive when moving into region II. 



comparable to ^. We have discussed the paranematic 
texture between two colloidal particles, showing that a 
defect of the paranematic director, in the form of a ring 
of topological charge —1/2 appears in the midplanc of 
the two particles. The diameter of the ring is approxi- 
mately equal to the distance between the centers of the 
two particles, which is easily explained by supposing that 
at lowest order the total paranematic tensor is the super- 
position of the paranematic tensors for the two isolated 
particles. The ring is surrounded by a sheath of biaxial 
order. 

Finally, using the Yukawa-like expression (|40|l for the 
paranematic interaction, we have examined the stabil- 
ity of a colloidal suspension by considering the interplay 
of the paranematic interaction with the standard DVLO 
interactions, i.e., van der Waals attractions and double- 
layer electrical repulsions. We have found that the stabil- 
ity of a colloidal suspension can be significantly affected 
by the paranematic interaction: the latter can either trig- 
ger flocculation or kinetically stabilize the suspension de- 
pending on the vicinity of the isotropic-nematic transi- 
tion of the nematogenic solvent and on the DLVO param- 
eters. These properties could be tested experimentally. 



APPENDIX A: DEFINITION OF THE 
SPHERICAL HARMONICS 



VII. CONCLUSIONS 

In this paper, following previous works [TsIITgL ITtL l20l 

|2]| , we have reconsidered the interactions between spher- 
ical colloids wetted by a corona of paranematic order in 
the isotropic phase of a nematogenic material. We have 
calculated the contribution to the total interaction that 
is mediated by the paranematic order, assuming a ra- 
dial anchoring of the director. Our results, based on a 
quadratic approximation of the free-energy density, are 
exact in the limit of weak induced paranematic order if 
the temperature is not too close to the nematic-isotropic 
phase transition. We have obtained analytical results in 
the asymptotic regime where the distance between the 
colloids is large with respect to the coherence length ^ 
of the nematic order, and numerical results for any sep- 
aration. For large separations, the interaction follows 
the Yukawa- like potential given by Eq. H40|l . which can 
be explained on simple grounds. Our numerical results 
show that, for colloids large compared to ^, this sim- 
ple potential is actually a very good approximation even 
up to separations of the order of f . We obtain attrac- 
tive interactions for colloids much larger than ^ and the 
possibility of short-range repulsions for colloids of size 



Among the various possible definitions of the spherical 
harmonics, we use the form 



where the P^"^{cosO) are associated Legendre functions 
of the first kind: 

Pr{t)^t^^r^tT'' !L_^r^t^Y , (A2) 

The Legendre functions obey the orthogonality relations 



P/"(i)P,-"(t) dt = , \ ZrJ' J u'^ (A3) 



{l + l/2){t-m) 



where is the Kronecker delta. These orthogonality 
relations imply the orthonormality of the spherical har- 



r^m(^, C^)Yl,^, {0, 0) sin d0d^ ^ Su'Sr, 



(A4) 
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